//
//  Simulation.cpp
//  phy_four
//
//  Created by Snir Gazit on 3/12/12.
//  Copyright (c) 2012 __MyCompanyName__. All rights reserved.
//

#define NMODEL 2.

#include <iostream>
#include "Simulation.h"
#include "boost/chrono.hpp"
#include <gsl/gsl_sf_hyperg.h>
#include <gsl/gsl_sf_gamma.h>

typedef boost::chrono::system_clock::time_point TimePoint;
typedef boost::chrono::duration<double, boost::ratio<3600>> Hour;
typedef boost::chrono::duration<double, boost::micro> MicroSec;

double calcW(double mu, double U, int m) {
	double c1, c2, c3, md = double(m), arg;
	arg = 0.25 * pow(mu, 2) / U;
	double f1 = 0.25 * (NMODEL + 2. * (1. + md));
	double f2 = 0.25 * (NMODEL + 2. * md);
	c1 = 0.25 * pow(U, -f1);
	c2 = mu * gsl_sf_gamma(f1) * gsl_sf_hyperg_1F1(f1, 3. / 2., arg);
	c3 = sqrt(U) * gsl_sf_gamma(f2) * gsl_sf_hyperg_1F1(f2, 1. / 2., arg);
	return c1 * (c2 + c3);
}
void calcWeight(double muIn, double UIn, double g, int maxBond,
		ostream &outFile, BoltzWeightType &weight) {
	double c, n, mu = muIn * g, U = UIn * g;
	c = calcW(mu, U, 0);
	int m = 1;
	for (auto w = weight.begin(); w != weight.end(); w++) {
		n = calcW(mu, U, m++);
		*w = (n / c);
		c = n;
	}
//	outFile << "BoltzRatiow: \n";
//	for (auto it = weight.begin(); it != weight.end(); it++)
//		outFile << *it << " ";
//	outFile << endl;
}
// Worm Alg
void Worm::init(QmcParameters &qmcParameters) {
	//Init Params
	qmcParameters_ = qmcParameters;
	// OutFile
	String fileName = qmcParameters_.outPath.getVal() + "/out.txt";
	outFile_.open(fileName.c_str());
	BoltzWeightType boltzRatio_;
	// Calc Weights
	if (qmcParameters_.simType.getVal() == "XY") {
		//boltzRatio_ = BoltzWeightType(qmcParameters_.maximalBonds.getVal(), 1.);
		boltzRatio_.fill(1.);
		qmcParameters_.g.setVal(1 / (2. * qmcParameters_.g.getVal()));
	} else {
		qmcParameters_.g.setVal(1. / (qmcParameters_.g.getVal()));
		calcWeight(qmcParameters_.mu.getVal(), qmcParameters_.U.getVal(),
				qmcParameters_.g.getVal(), qmcParameters_.maximalBonds.getVal(),
				outFile_, boltzRatio_);
	}
	// Init Lattice
	lat_.init(qmcParameters_.Llist.getVal());
	lat_.getIra()->totalBonds_++;
	// Init Rng
	eng_.seed(qmcParameters_.seed.getVal());
	// Init Moves
	shiftMove_.init(&eng_, &lat_, qmcParameters_, boltzRatio_);
	jumpMove_.init(&eng_, &lat_, qmcParameters_, boltzRatio_);
	// Measurements
	measures_.init(&lat_, qmcParameters_, &boltzRatio_, &outFile_);
	// Open Close
	UniIntParams param(0, 1);
	openClose_.init(&eng_, param);
	// Read From File
	if (qmcParameters_.initFromFile.getVal())
		this->read();
	int totalNumberOfMeasure = (int) pow(
			(double) qmcParameters.binSize.getVal(),
			qmcParameters.numOfBins.getVal());
	measures_.updateMinNumOfMeasure(totalNumberOfMeasure);
}
Worm::~Worm() {
	outFile_.close();
}
void Worm::nextMove() {
	// Check in Closed
	if (lat_.isClosed()) {
		// Jump Or Shift
		int openClosed = openClose_.draw();
		if (openClosed == 0)
			jumpMove_.makeMove();
		else
			shiftMove_.makeMove();
		// Shift
	} else {
		shiftMove_.makeMove();
	}
}

void Worm::thermalize() {
	int numOfThermal = qmcParameters_.numOfThermal.getVal();
	for (int i = 0; i < numOfThermal; i++) {
		this->nextMove();
	}
	int numOfClosed = 0;
	for (int i = 0; i < numOfThermal; i++) {
		this->nextMove();
		if (lat_.isClosed())
			numOfClosed++;
	}
	double pOpen = double(numOfClosed) / double(numOfThermal);
	int closedRate = int(
			double(qmcParameters_.rateFactor.getVal())
					* double(lat_.getNumOfSites()) * pOpen);
	measures_.updateClosedRate(closedRate);
	outFile_ << "Closed Rate - " << closedRate << endl;
}

bool Worm::run() {
	bool endTime, dumpFlag;
	while (measures_.checkNotDone()) {
		this->nextMove();
		measures_.DoMeasurement(endTime, dumpFlag);
		if (dumpFlag) {
			this->dump();
		}
		if (endTime) {
			this->dump();
			break;
		}
	}
	return !endTime;
}

void Worm::dump() {
	String outPath = qmcParameters_.outPath.getVal();
	measures_.dump(outPath);
	lat_.dump(outPath);
	eng_.dump(outPath);
}
void Worm::read() {
	String outPath = qmcParameters_.outPath.getVal();
	measures_.read(outPath);
	lat_.read(outPath);
	eng_.read(outPath);
}

// WormLoop Algorithm

// Worm Alg
void WormLoop::init(QmcParameters &qmcParameters) {
	//Init Params
	qmcParameters_ = qmcParameters;
	// OutFile
	string runNum = qmcParameters.runNum.getVal();
	// OutFile
	String fileName = qmcParameters_.outPath.getVal() + "/out" + runNum
			+ ".txt";
	outFile_.open(fileName.c_str());
	outFile_ << "Stating N=2 Simulation - Run Number - " + runNum << endl;
	if (qmcParameters_.simType.getVal() == "XY") {
		//weight = BoltzWeightType(qmcParameters_.maximalBonds.getVal(), 1.);
		boltzRatio_.fill(1.);
		qmcParameters_.g.setVal(1 / (2. * qmcParameters_.g.getVal()));
	} else {
		qmcParameters_.g.setVal(1. / (qmcParameters_.g.getVal()));
		calcWeight(qmcParameters_.mu.getVal(), qmcParameters_.U.getVal(),
				qmcParameters_.g.getVal(), qmcParameters_.maximalBonds.getVal(),
				outFile_, boltzRatio_);
	}

	// Init Lattice
	lat_.init(qmcParameters_.Llist.getVal());
	// Init Rng
	eng_.seed(qmcParameters_.seed.getVal());
	// Init Moves
	CreateLoopMove_.init(&eng_, &lat_, qmcParameters_, boltzRatio_);
	// Measurements
	measures_.init(&lat_, qmcParameters_, &boltzRatio_, &outFile_);
	// Read From File
	if (qmcParameters_.initFromFile.getVal())
		this->read();
	int totalNumberOfMeasure = (int) pow(
			(double) qmcParameters.binSize.getVal(),
			qmcParameters.numOfBins.getVal());
	measures_.updateMinNumOfMeasure(totalNumberOfMeasure);
}
WormLoop::~WormLoop() {
	outFile_.close();
}
void WormLoop::simulatedAnnealing() {
	outFile_ << " Startting Simulated Anneling " << endl;
	double gHighT = 1. / qmcParameters_.gHighT.getVal();
	int numOfSMA = qmcParameters_.numOfGSA.getVal();
	double dg = (qmcParameters_.g.getVal() - gHighT) / double(numOfSMA - 1);
	acc wormLength;
	double nLength;
	int numOfThermal = qmcParameters_.numOfThermal.getVal();
	vector<double> gSA;
	for (int i = 0; i < qmcParameters_.numOfGSA.getVal(); i++)
		gSA.push_back(gHighT + dg * double(i));
	for (auto it_g = gSA.begin(); it_g != gSA.end(); it_g++) {
		double cg = *it_g;
		setCoupeling(cg);
		outFile_ << "Current g -" << cg << endl;
		wormLength = acc();
		for (int i = 0; i < numOfThermal; i++)
			wormLength(double(CreateLoopMove_.makeMove()));
		nLength = ba::mean(wormLength);
		numOfThermal = int(
				100 * double(lat_.getNumOfSites()) / double(nLength));
		outFile_ << "Worm Length - " << nLength << " numOfThermal "
				<< numOfThermal << endl;
	}
	numOfThermal = int(1000 * double(lat_.getNumOfSites()) / double(nLength));
	wormLength = acc();
	for (int i = 0; i < numOfThermal; i++)
		wormLength(double(CreateLoopMove_.makeMove()));
	nLength = ba::mean(wormLength);
	int closedRate = int(
			double(qmcParameters_.rateFactor.getVal())
					* double(lat_.getNumOfSites()) / nLength);
	measures_.updateClosedRate(closedRate);
	outFile_ << "Closed Rate - " << closedRate << endl;
}

void WormLoop::setCoupeling(double g) {
	if (qmcParameters_.simType.getVal() == "HEISENBERG")
		boltzRatio_.fill(1.);
	else
		calcWeight(qmcParameters_.mu.getVal(), qmcParameters_.U.getVal(), g,
				qmcParameters_.maximalBonds.getVal(), outFile_, boltzRatio_);
	CreateLoopMove_.setCoupeling(boltzRatio_, g);
}
void WormLoop::thermalize() {
	int numOfThermal = qmcParameters_.numOfThermal.getVal();
	namespace ba = boost::accumulators;
	typedef ba::accumulator_set<double, ba::stats<ba::tag::mean> > acc;
	acc wormLengthMA;
	for (int i = 0; i < 10 * numOfThermal; i++) {
		wormLengthMA(CreateLoopMove_.makeMove());
	}
	double cLength = ba::mean(wormLengthMA), nLength;
	bool conv;
	do {
		wormLengthMA = acc();
		for (int i = 0; i < numOfThermal; i++) {
			wormLengthMA(double(CreateLoopMove_.makeMove()));
		}
		nLength = ba::mean(wormLengthMA);
		outFile_ << "Worm Length - " << nLength << endl;
		double ratio = fabs(1. - (cLength / nLength));
		conv = ratio < 0.01;
		cLength = nLength;
		//	wormLengthMA = acc();
	} while (!conv);
	//CreateLoopMove_.setAvWormLength(2 * ba::mean(wormLength));
	outFile_ << "Final Worm Length - " << cLength << endl;
	int closedRate = int(
			double(qmcParameters_.rateFactor.getVal())
					* double(lat_.getNumOfSites()) / cLength);
//	closedRate = 1;
	measures_.updateClosedRate(closedRate);
	outFile_ << "Closed Rate - " << closedRate << endl;
}

bool WormLoop::run() {
	bool endTime, dumpFlag;
	namespace ba = boost::accumulators;
	ba::accumulator_set<double, ba::stats<ba::tag::mean> > wormLength;
	while (measures_.checkNotDone()) {
		wormLength(double(CreateLoopMove_.makeMove()));
		measures_.DoMeasurement(endTime, dumpFlag);
		if (dumpFlag) {
			this->dump();
		}
		if (endTime) {
			this->dump();
			break;
		}
	}
	outFile_ << " Worm Length - " << ba::mean(wormLength) << endl;
	outFile_ << " Erase Prob - " << CreateLoopMove_.getEraseProb() << endl;
	return !endTime;
}

void WormLoop::dump() {
	String outPath = qmcParameters_.outPath.getVal();
	measures_.dump(outPath);
	lat_.dump(outPath);
	eng_.dump(outPath);
	String fileName = outPath + "/sim_params.xml";
	// Write to XML
	pugi::xml_document doc;
	pugi::xml_node rootNode = doc.append_child("SIM_PARAMS");
	// Parameters
	pugi::xml_node node = rootNode.append_child("CLOSED_RATE");
	node.append_child(pugi::node_pcdata).set_value(
			(boost::lexical_cast<string>(measures_.getClosedRate())).c_str());
	doc.save_file(fileName.c_str());
//	doc.save(std::cout);
//	cout << endl;
}

void WormLoop::read() {
	outFile_ << "Starting Reading From File" << endl;
	String outPath = qmcParameters_.outPath.getVal();
	measures_.read(outPath);
	outFile_ << "Measures Loaded" << endl;
	lat_.read(outPath);
	int totalBonds;
	try {
		bool testConf = lat_.testConfiguration(totalBonds);
		outFile_ << "Total Bonds Calculated - " << totalBonds
				<< " total Bonds previous sim - " << lat_.getTotalBonds()
				<< endl;
		if (totalBonds != lat_.getTotalBonds())
			throw "Total Bonds after Read don't match";
		if (!testConf)
			throw "Bad Load of Configuration";
	} catch (string &e) {
		outFile_ << e;
		exit(1);
	}
	outFile_ << "Lattice Loaded" << endl;
	eng_.read(outPath);
	String fileName = outPath + "/sim_params.xml";
	// Read From XML
	pugi::xml_document doc;
	try {
		pugi::xml_parse_result result = doc.load_file(fileName.c_str());
		if (!result)
			throw result;
	} catch (pugi::xml_parse_result &result) {
		std::cout << "Error description: " << result.description() << "\n";
		exit(1);
	}
	pugi::xml_node rootNode = doc.child("SIM_PARAMS");
	measures_.updateClosedRate(readXmlData<int>(rootNode, "CLOSED_RATE"));

}
